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Image Demosaicing and Enhancement System 
Field of the Invention 

The present invention relates to digital cameras, and more particularly, to an improved 
method for converting data from a camera sensor to a color image. 

Background of the Invention 

A digital color image usually consists of an array of pixel values representing the intensity of 
the image at each point on a regular grid. Typically, three colors are used to generate the 
image. At each point on the grid the intensity of each of these colors is specified, thereby 
specifying both the intensity and color of the image at that grid point. 

Conventional color photography records the relevant image data by utilizing three 
overlapping color sensing layers having sensitivities in different regions of the spectrum 
(usually red, green, and blue). Digital cameras, in contrast, typically utilize one array of 
sensors in a single "layer". 

When only one sensor array is used to detect color images, only one color may be detected at 
any given sensor location. As a result, these sensors do not produce a color image in the 
traditional sense, but rather a collection of individual color samples, which depend upon the 
assignment of color filters to individual sensors. This assignment is referred to as the color 
filter array (CFA) or the color mosaic pattern. To produce a true color image, with a full set 
of color samples (usually red, green and blue) at each sampling location, a substantial amount 
of computation is required to estimate the missing information, since only a single color was 
originally sensed at each location in the array. This operation is typically referred to as 
"demosaicing". 
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To generate the missing information, information from neighboring pixels in the image 
sensor must be used. A number of algorithms have been put forward in an attempt to provide 
the missing information while minimizing artifacts resulting from the estimation process. 
The simplest algorithms interpolate the sensor data from like color sensors to provide the 
missing information. These algorithms treat the red sensors as being independent from the 
green sensors, and so on. To provide a red value at a given location, the values measured by 
the red sensors in the region of that location are interpolated. This approach requires that the 
image be low-pass filtered. Such filtering reduces the image resolution below the pixel 
resolution of the underlying sensor array. This lost resolution cannot be recovered. 

To avoid this loss in resolution, less aggressive optical low-pass filtering is used in some 
higher-end cameras. However, in such systems, the color sensors may no longer be treated as 
independent. For example, Wober, et al. (U.S. Patent 5,475,769) describe a method for 
generating the missing color information by computing a weighted average of the pixel 
values in the neighborhood of the pixel whose missing color information is being computed. 
This method weights values from all of the color sensors, not just the color being 
reconstructed. However, even this approach leaves much to be desired since it utilizes one set 
of weights for all images. 

A single pixel array may be viewed as consisting of a number of separate planes of pixels in 
which each plane has sensors for the same color. Since the pixels do not overlay, the sensors 
in the various planes are at different locations. Systems that take weighted averages across 
more than one plane make use of the statistical dependencies between these sample locations. 
In effect, the blurring of an image by the camera optics allows an image edge that falls on one 
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color plane precisely on the sensors of that plane to also be seen in the other color planes 
because the image is spread by blurring onto the sensors in the other color plane. Since the 
statistical dependencies between the various color planes depend on the amount of blur 
introduced by the camera optics, an optimal algorithm must take into account the physical 
camera settings. Accordingly, a single set of weight functions will not provide an optimal 
estimation of the missing information. 

The statistical dependencies also depend on the source of illumination. Different illumination 
sources have different spectra. The pixel filters have broad pass-bands centered at the red, 
green, and blue wavelengths. In the absence of any image blurring, the response of any given 
pixel is determined by its color filter, the reflectivity of the corresponding point in the scene 
being photographed, and the light spectrum incident on that point from the illumination 
source. The blurring provided by the camera optics mixes the light between the pixels. 
Hence, the statistical dependencies, in general, depend both on the illumination source and 
the camera optics. Prior art methods for converting the pixel array data to a fully sampled 
color digital image do not take the illumination source into account. 

Broadly, it is the object of the present invention to provide an improved image processing 
method for converting data from a pixel array having non-overlapping sensors to a fully 
sampled digital image. 

It is a further object of the present invention to provide a conversion method that corrects for 
the camera's optical system. 



It is a still further object of the present invention to provide a conversion method that corrects 
for the source of illumination. 

These and other objects of the present invention will become apparent to those skilled in the 
art from the following detailed description of the invention and the accompanying drawings. 

(K-> Summary of the ■ Invention ' 

The present invention is a method for operating a data processing system to generate a second 
image from a xirst image. The first image includes a two dimensional array of pixel values, 
each pixel value (Wresponding to the light intensity in one of a plurality of spectral bands at a 
location in the first miage. The method utilizes a linear transformation of a vector derived 
from super input pixels to obtain a vector that includes at least one super output pixel. The 
super input pixels are defined by separating the pixels of the first image into a plurality of 
input image planes. Each input image plane has an identical number of pixels within a 
normalized horizontal and vertical sampling interval as the other input image planes. All 
pixels in a given input image plane correspond to tt^same spectral band as the other pixels in 
that input image plane. Each super input pixel is a ve^tor^of dimension P, where P is the 
number of the input image planes, each exponent of that vector being an input pixel from a 
corresponding input image plane. Similarly^, set of output image planes is defined, each 
pixel in a given output image plane representin^the intensity of the second image in one of a 
plurality of spectral bands at a corresponding point in the second image. Each super output 
pixel is a vector of dimension Q, where Q is the number of the output image planes, each 
component of that vector being a pixel from a corresponding output image plane. In the 
preferred embodiment of the present invention, the linear transformation depends on the 



properties of\the optical system and the illumination source used to generate the first image. 
The linear transformation can also be varied to take into account the contents of the scene 
captured in the firs^mage and the desired output format of the second image. 

Brief Description of the Drawing s 

Figure 1 illustrates the separation of an image taken with an image sensor having a repeating 
2x2 pattern into image planes according to the present invention. 

Figure 2 illustrates a portion of a sensor array and the input pixels in the sensor array which 
contribute to a particular intermediate input vector. 

Figure 3 illustrates a portion of an output RGB (red, green, blue) image and the pixels in the 
RGB output image that correspond to the intermediate output vector shown in Figure 2. 

Detailed Description of the Invention 

The method of the present invention may be applied to any color-sampling device that 
acquires its samples from a sensor array that can be decomposed into a plurality of image 
planes, which satisfy two conditions. First, each image plane must have an identical number 
of samples within a normalized horizontal and vertical sampling interval; however, the 
various image planes may be arbitrarily displaced relative to one another. Second, all 
samples in a given image plane must have identical color properties; however, multiple image 
planes can have the same color properties. 



+ 



3 S 



fii 

D 

i y 

I s § 



10 



15 



These conditions are satisfied by any image sensor having a sensor pattern that is constructed 
)y repeating a kernel of sensing elements. For example, one common image sensor array is 
based on the IJayer pattern, which generated by repeating a 2x2 sensor array kernel having 
two green sensors, one red sensor, and one blue sensor. This pattern is shown in Figure 1 at 
10. The kernel isshown at 12. Such an image sensor may be viewed as having 4 planes 
shown at 14-17, two green planes 14 and 17, one red plane 16, and one blue plane 15. The 
sampling interval is the area originally occupied by the kernel. Each of the planes is offset 
with respect to the other planes. It can be shown that any regular sampling lattice can be 
decomposed into a set of imaget planes satisfying the above conditions. 

To simplify the following discussion, vector notation will be utilized. Vectors and matrices 
will be shown in bold print to distinguish them from scalar quantities. The measured 
intensity values in each image plane will be denoted by Xp[n l5 n 2 ]. Here, n x and n 2 are indicies 
which denote the position of the pixel in the p^ image plane and Xp is the intensity value 
measured for that pixel. The quantity [n u n 2 ] is a two-dimensional integer valued vector 
which will be denoted by n. The entire set of image planes can then be represented as a set of 
vectors x[n] where 



x[n] = 



x,[/i,,n 2 ]" 
x 2 [n,,« 2 ] 

L x,[»,,« 2 ]. 



(1) 
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The output image can likewise be represented as a set of vectors defined on a different set of 
image planes. Typically, the goal of the demosaicing algorithm is to generate a set of 
regularly spaced pixels in three color planes (red, green, and blue). Denote the intensity in 



the i* color plane by y^n,,^]. Then the output pixels can also be represented by a set of 
vectors. 



y[m] = 



y 2 [m l9 m 2 ] 
y Q [m l9 m 2 ] 



(2) 



In the demosaicing case, Q is typically 3; however, different Q values may be utilized. For 
example, an image that is to be printed on a color printer utilizing 4 dyes could be generated 
directly by the method of the present invention utilizing a representation in which Q=4. 

In general, the output image will have a spatial resolution that is different from the input 
image. The input image may be viewed as consisting of a set of "super pixels", x[n]. 
Likewise, the output image is a set of pixels y[m]. The number of output pixels in the 
vertical and horizontal directions corresponding to each input pixel will be denoted by A, and 
X 2 , respectively. In the case of the Bayer pattern discussed above, the demosaicing task is 
usually understood as having A,, = X 2 = 2. That is, one attempts to construct one output (RGB) 
pixel for each physical sensor in the input array. 

In the method of the present invention, the output pixels are related to the input pixels by a 
linear operator that operates on vectors derived from x[n] and y[m]. These intermediate 
vectors take into account the difference in resolution and the fact that each output pixel 
depends on more than one input super pixel. The intermediate vector corresponding to y[m] 
will be denoted by <^[n] and has the same sampling density as x[n]: 



y[A(n + 8,)] 
y[A(n + 8 2 )] 

_y[A(n + 8 XiX2 )]_ 



Here, the matrix A is defined by 



A = 



\ 0 
0 X, 



(3) 



(4) 



In the case of the Bayer pattern, 8,=[0,0], 8,=[ 1/2,0], 8,=[0,l/2], and 8,=[l/2,l/2]. The 
vectors ^[n] will be referred to as the output polyphase components in the following 
discussion. 
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In tne method of the present invention, it is assumed that each polyphase output vector 
depends\n a finite number of input super pixels. In general, the input super pixels that 
contribute tcNa particular polyphase output vector £[n] will be located in a neighborhood 
around [n]. As win be explained in more detail below, the precise pixels will depend on the 



nature of the camera\and 



imaging optics. The input super pixels that contribute to the 
polyphase output vector at\n may be identified by a set of displacement vectors k„ k 2 , . . ., k K . 
That is, <^[n] depends on xfn+kj, x[n+k 2 ], x[n+k K ]. In the method of the present 
invention, £[n] is assumed to be linearly dependent on the input super pixels. In the preferred 
embodiment of the present invention, the set of displacement vectors k„ k 2 , k K is 
independent of n, and is arranged in a k x x\k 2 rectangular grid 
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The linear relationship can be most easily defined in terms of a vector t\n] which includes all 
of the super pixels on which the output polyphase vector ^[n] depends, i.e., 



"x[n + k,] 
x[n + k K ]_ 



(5) 



In terms of this vector, the relationship between C[n] and £[n] may be written as matrix 
multiplication: 

C[n]=T^[n] (6) 



where T is a (QX 1 X 2 )x(PK 1 K 2 ) matrix. 



Refer now to Figures 2 and 3 which illustrate the relationships between the output pixels 
Yj[N,M] ? the input pixels XjflNfjM] and the two intermediate vectors defined above for the 
Bayer sensor pattern. Figure 2 illustrates a portion of a sensor array and the input pixels in 
the sensor array which contribute to ^[N,M]. Figure 3 illustrates a portion of an output RGB 
image and the pixels in the RGB output image that correspond to ^[N,M] and which are 
computed from the pixels shown in Figure 2 by the matrix multiplication operation shown in 
Eq. (6). 



The matrix, T, depends on a number of factors. Some of these are fixed for a particular 
imaging device and some depend on the particular manner in which the imaging device is 
being utilized. For example, the physical properties of the sensing array such as the spectral 
sensitivity of the pixels, the mosaic pattern, and the number of pixels typically do not vary 
from image to image. In contrast, the optical properties of the imaging device such as the 
lens settings on a camera (f number and zoom) may vary from image to image. In addition, 
the spectral properties of the illumination source may vary from image to image (daylight, 
flash, incandescent light, etc.). 
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In addition, the statistics of the image being captured may be taken into account through T. 
For example, images or portions of images having a high content of vertical and horizontal 
edges can be processed with a different matrix than images lacking such features and thereby 
improve the output image quality. 

In cameras that have a variable output format, the resolution of the final picture can be set 
using a different T matrix. Alternatively, a single T matrix may be utilized for all resolutions 
and then the desired output image determined by re-sampling the fixed resolution image. 
Similarly, the number of output color planes may be altered by using different T matrices or 
by resampling a single color format to generate an alternate color representation. In general, 
properties that alter the dimension of the T matrix are preferably handled by using a fixed T 
matrix and then re-sampling the final image in low cost imaging devices. 

If the number of different T matrices is relatively small, the coefficients of the T matrix can 
be determined by training the system on known images. For each possible T matrix, images 
of a number of known scenes are taken using the imaging device. The coefficients of the T 
matrix are then computed so as to minimize the difference between the image computed from 
the sensor input and the known scene images. Such optimization computations are well 
known to those skilled in the art, and hence, will not be discussed in detail here. 

If the variation in some parameter such as f-number is relatively smooth, the T matrices need 
only be computed for a discrete number of values of the parameter. The correct T matrix for 
the non-computed variable parameter values can then be computed by interpolation of the 
computed T matrices. 
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Model-Based Computation of T 

As noted previously, in some circumstances it may be possible to compute appropriate 
matrices, T, from training images. Unfortunately, this approach is limited to appli- 
cations in which the number of different imaging conditions and hence the number 
of different T matrices which could be required is relatively small. The purpose of 
the material presented here is to describe a method for directly computing T for an 
arbitrary imaging device (i.e. arbitrary color sensitivities, sensor locations and optical 
characteristics) and under arbitrary illumination, subject to a particular statistical 
model for the underlying image, which has been found to give particularly good re- 
constructed image quality. As will be seen, the statistical image model is governed 
by only a few parameters. In more advanced applications, these parameters may be 
adjusted, either locally, or globally, to match statistical properties of the image, such 
as edge orientation, which can be estimated by various methods. 

Image Formation Model 

This section describes the parameters of the image formation process which maps the 
original scene into the source image super-pixels, x[n]. The image formation model 
depends upon deterministic quantities which can, at least in theory, be measured. 
These quantities are 

• The scene illuminant spectral power density, /(A). 

• The color spectral response functions, r p (A), for each input image plane, p. 

• The Point Spread Function, h p (X, s), associated with the combined effects of the 
optical transfer function and sensor integration behaviour for input image plane 
p. Here s = [^i,^] is the spatially continuous argument of the Point Spread 
Function (PSF), at each wavelength, A. Note that the PSF's implicitly include 
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the effect of relative displacements between the different input image planes. 
In the following, the PSF is referenced only through its Fourier Transform, 
/tp(Acj), where the spatial frequency vector, u = [i*>i,u>2]> is normalized so that 
= u? 2 = 7r at the Nyquist frequency of the input super-pixels. Thus, the 
Nyquist frequency of the original sensor array corresponds tou = [Ai7r, X 2 tt]. 

Rather than modeling the image formation process directly in terms of the desired out- 
put image super-pixels, it is helpful to choose an intermediate representation in terms 
of surface spectral-reflectance, since this is well known to be better behaved from a 
statistical perspective than the scene radiance itself and hence more amenable to the 
statistical modelling described in the next section. Specifically it is helpful to assume 
that the spectral reflectance of the original scene can be perfectly represented as a lin- 
ear combination of a limited number of fixed basis functions, 61(A), 62(A), . . . , 65(A), 
where S is usually chosen to be three or four, but may be larger if desired. The 
actual output vectors, y[m], may be expressed in terms of the intermediate spectral 
reflect ance vectors, z[m] as 

y[m] = T out • z[m] 

f JS° d x (\)l(\)bi{\)d\ JT d x (X)l(X)b 2 (X)dX ... d 1 (X)l(X)b s (X)dX } 
JS° d 2 (X)l(X)b x (X)dX f£° d 2 (X)l(X)b 2 (X)dX ... f£° d 2 (X)l(X)b s (X)dX 

{ fo° d Q (X)l(X)b 1 (X)dX f™ d Q (X)l(X)b 2 (X)dX ... d Q (X)l(X)b s (X)dX ) 

where d q (X) is the spectral response of the q'th display spectral response function. For 
example, if the objective is to recover an XYZ image, then Q should be set to 3 and 
dx(X) through d 3 (X) should be set to the standard 1931 CIE tri-stimulus functions. 
As another example, if the objective is to recover an image with the same color 
characteristics as the different color filters on the physical sensor array, then Q should 
be set to the number of unique input response functions, r p (A), and there should 
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be a one-to-one correspondence between these unique r p (A) and the d q (\). In this 
framework, the chief goal is to compute the (SXi A 2 ) x (PKiK 2 ) reconstruction matrix, 
T re / which maps the neighbourhood of input super-pixels, £[n], to the corresponding 
spectral-reflectance super-pixel, 

( z[A(n + *,)] ^ 

z[A(n + £ 2 )] 
C [n] = . 

< z[A(n + 6 Xl x 2 )] j 

The final (QX1X2) x (PK1K2) reconstruction matrix is then formed by simple matrix 
multiplication: 

' T ont 0 0 ••• 0 0 ^ 
0 T out 0 • • 0 0 



0 



0 



0 T out J 



ref 



(7) 



The linear image formation model may now be expressed compactly as 

x(w) = u(w) + £ HMziA- 1 ^) 

where 



• x(u>) is the Discrete Space Fourier Transform of the input image, x[n]; 

• is the Discrete Space Fourier Transform of the sampling noise vector se- 
quence; 

• z(w) is the Discrete Space Fourier Transform of the spectral reflectance vector, 
z[m]; 



• 
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f2 a (o;), is the set containing all AiA 2 aliasing frequencies associated with the 
sampling of the high resolution output image on grid [m] onto the input super- 
pixel grid [n], for each lj € [— 7r,7r] 2 ; and 



H(u>) is the P x S image formation matrix, 



H(u) = 



( / 0 oo r 1 (A)/(A)6 1 (A)^(A,u;)dA 
/ 0 °°r 2 (A)/(A)6 1 (A)^ 2 (A,a;)rfA 



/ 0 °° r P (\)\{\)b x {\)hp{\^)d\ 



/ 0 oo r 1 (A)/(A)6 s (A)A 1 (A,u;)dA 
/ 0 oo r 2 (A)/(A)6 s (A)^(A,u;) t ZA 



Statistical Model 



In order to compute an appropriate solution to the image reconstruction problem, 
it is necessary to introduce a statistical model for the sampling noise, and the 

spectral reflectance, z(u>). In this discussion, Wide Sense Stationary Gaussian models 
are assumed, which are characterized entirely by the covariance matrices, 

C„(w) = E \i>{u) ■ u{uY] 

and 

C,(w) = E [z(w) • z(w)*] 

The noise covariance matrix will usually be a constant, C„(w) = cr 2 I, for all u>, 
corresponding to white noise, but other models may be used if appropriate. 

The following parametric model is used for the reflectance covariance matrix, 

C,(u;) = C B 0 .|ra;|"* (8) 

where C z ° is a constant 5x5 covariance matrix, p is a frequency roll-off parameter, 
which is usually selected in the range 20 to 30 dB/decade, and T is a 2 x 2 "shape 
matrix". The terms in the above expression which follow the constant covariance 
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matrix, C z °, describe a scalar envelope function whose contours are elipses in the 
frequency domain. The orientation and aspect ratio of these eliptic contours may 
be explicitly controlled by means of the T matrix. For a circular cross-section, the 
identity matrix, T = I, may be used. 

The statistical model represented by C z plays an extremely important role in deter- 
mining the quality of the final reconstructed images. The parametric model described 
above may be justified on a number of grounds; most notably, the model is scale- 
invariant, which means that on average the statistics of scenes should not depend 
upon how far the camera is located from the scene. This scale- invar iance property is 
important because in practical imaging applications, information about the absolute 
scale of objects in the scene is rarely available. Also, there is significant empirical 
evidence for this scale invariance property in natural scenes, with a frequency roll-off 
factor, p, of about 20 dB/decade. 

Efficient Computation of T 

As mentioned above, the key objective is to compute the (S\i\ 2 ) x (PKiK 2 ) matrix, 
T re £, from which T is easily recovered via equation (7). The ensuing discussion 
concerns the derivation of an optimal Linear Minimum Mean Squared Error (LMMSE) 
estimator, T re f, subject to the models described in the previous two sections. The 
formula for such an estimator is well-known. Specifically, 

T re f = Z • X- 1 (9) 
where Z is the (SX1X2) x (PKiK?) cross-covariance matrix, 

Z = E [C'[n] - £[n] 4 ] 
and X is the (PKiK 2 ) x (PK X K 2 ) dimensional auto-covariance matrix, 

X = E [tfn] • £[n]«] 
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In fact, X has the following Toeplitz block-Toeplitz structure, 



X = 



X[0] 

x[-i] 



x[i] • •• x[k 2 -i] N 



X[0] 



V X[l-if 2 ] X[2-/iT 2 ] 
where each block, X[/ 2 ], has the Toeplitz form 

( xjo,z 2 ] xji,/ 2 ] 
x[-i,/ 2 ] xjo,/ 2 ] 



X[/r 2 - 2] 



X[0] ) 



X[#i-U 2 ] 
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\ Xjl - , / 2 ] X[2 - #i , h] ••■ X[0, / 2 ] / 

and each sub-block, Xj/ 1? / 2 ], is a P x P source super-pixel covariance matrix, given 
by 

= E [x[n] • x[n + 1]'] 

The (S'A 1 A 2 ) x (PKiK 2 ) matrix, Z, also has a doubly-nested block structure. Specif- 
ically, 

where 

■w-(i[-i^Hi[ i -iTj^]-*m^]) 

and the (SA1A2) x P sub-block matrices, £[1], are given by 

1(1] = E [C[n] • £[n + 1]'] . 

In order to compute T re j, then, it is sufficient to compute the matrices for 
[1] = ['i,^] in the range -I<i < U < I<i and the matrices £[1], for [1] = [l u l 2 ] in the 
range — [^j < h < [^2^] > after which the contents of X and Z may be filled in 
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and used to evaluate equation (9). The key to efficient computation of T re f lies in 
efficiently computing the matrices, and 

It turns out that these matrices may be efficiently computed by exploiting Parseval's 
Relationship. Specifically, 



yZir) J-it J— TV 



and 



= [_l] = (2^ I-**" 1 /' d ^ cVIC C*( w ) 



(10) 



(11) 



where the frequency auto- and cross-covariance matrices, C x (u;) and C£/ X (u;) are 
found from 

C x («) = C„(w) + £ (HK)C z (u, 0 )H(u, a )*) 

w 0 6n 0 (u;) 

and 

C<.x(«) = £ (*("a)C z (a, 0 )H(u, a r) 

Here, <&(u>) is the (5AiA 2 ) x S matrix of phase shifts, 

/ e i^ £ 5il(5) \ 



(12) 



(13) 



*(u;) = 



e i^ t ^A 1 A 2 i(5) j 

corresponding to the relative dispacements of each of the output polyphase compo- 
nents. 

In order to compute the matrices, and 2L[1], C^ x (uj) and C Cx (u>) are evaluated at 
a finite number of frequencies, u> G [— 7r,7r] 2 and then the Inverse Fourier Transform 
(IFT) integrals of equations (10) and (11) are approximated numerically. There are 
various approaches to determining the best set of frequency points at which to evaluate 
X[l] and Z[l] and interpolating between these points during the numerical integration 
procedure, but these are beyond the scope of this brief discussion. 
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Various modifications to the present invention will become apparent to those skilled 
in the art from the foregoing description and accompanying drawings. Accordingly, the 
present invention is to be limited solely by the scope of the following claims. 



